% function[dx,dz] = eej_dx(w,h,x),
% dx = atan((2*w*h)/(h^2+x^2-w^2));
% dz = log10(((x+w)^2 + h^2)/((x-w^2 + h^2)))/(2*w);

CD = load('e:\manoj\projects\eej-induction\new_CD_prof.txt');
data = load('e:\manoj\projects\eej-induction\GROUND_5.TXT');

mu = 4*pi*10^-7;
h = 100*1e3;
%I_ion = 1e6;
I_ion = 65*1e3;
pp = 1;
prd = 2*3600;
% data1 = load('c:\manoj\projects\eej-induction\quebec.txt');
% data1(:,1) = data1(:,1)*1e3;
% data1(:,2) = 1./data1(:,2);
data1 = load('e:\manoj\projects\eej-induction\ak\thick_sigma_corrected.txt');
Y = cres(1,5,data1(:,1),data1(:,2),prd);
wm = sqrt(-1)*mu * 2.0*pi/prd;
p = conj(-1.0/(wm*Y(1)));

% P_fun = c_res
mu = 4*pi*10^-7;
for i = 1:length(data),
    dummy_e = 0;
    dummy_i = 0;
    dummy_b_e = 0;
    dummy_b_i = 0;
    for j = 1:length(CD),
                x = deg2km(distance(CD(j,1),78,data(i,1),78))*1e3;
                bin_size = 55.5419*1e3;
                h = 108*1e3;
                dummy_e = dummy_e + (mu*bin_size*CD(j,2)/(2*pi))*((h/(h^2+x^2)));
                dummy_i = dummy_i + (mu*bin_size*CD(j,2)/(2*pi))*((h+2*p)/((h+2*p)^2+x^2));
    end;
    Bxe(i) = dummy_e;
    Bxi(i) = dummy_i;
end;
                
%Ran the code on 9dec06.. The external and internal part of the EEJ - H mag
%profile at ground seems OK..



% pp=1;
% for x = -1000*1e3:10*1e3:1000*1e3,
% b_x_e(pp) = (mu*I_ion/(2*pi))*((h/(h^2+x^2)));
% b_x_i(pp) = (mu*I_ion/(2*pi))*((h+2*p)/((h+2*p)^2+x^2));
% pp = pp+1;
% end;

hold on;
plot(CD(:,1),real(Bxe)*1e9,'r');
hold on;
plot(CD(:,1),real(Bxi)*1e9,'k');